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Abstract. Progress in computing the spectrum of excited baryons and mesons in lattice QCD is de- 
scribed. Large sets of spatially-extended hadron operators are used. A new method of stochastically 
estimating the low-lying effects of quark propagation is utilized which allows reliable determina- 
tions of temporal correlations of both single-hadron and multi-hadron operators. The method is 
tested on the Tj . (7, o mesons. 
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We are currently carrying out computations of the excitation spectrum of QCD in 
finite volume with ab initio Markov-chain Monte Carlo path integrations on anisotropic 
space-time lattices. Our first results using two flavors of dynamical quarks were reported 
in Ref. [1], and our most recent results can be found in Ref. [2]. Such calculations are 
very challenging. Computational limitations cause simulations to be done with quark 
masses that are unphysically large, leading to pion masses that are especially heavier 
than observed. The use of carefully designed quantum field operators is crucial for 
accurate determinations of low-lying energies. To study a particular state of interest, 
the energies of all states lying below that state must first be extracted, and as the 
pion gets lighter in lattice QCD simulations, more and more multi-hadron states lie 
below the masses of the excited resonances. The evaluation of correlations involving 
multi-hadron operators contains new challenges since not only must initial to final time 
quark propagation be included, but also final to final time quark propagation must be 
incorporated. 

The use of operators whose correlation functions C{t) attain their asymptotic form as 
quickly as possible is crucial for reliably extracting excited hadron masses. An important 
ingredient in constructing such hadron operators is the use of smeared fields. Operators 
constructed from smeared fields have dramatically reduced mixings with the high fre- 
quency modes of the theory. Both link-smearing[3] and quark-field smearing[4] must 
be applied. Since excited hadrons are expected to be large objects, the use of spatially 
extended operators is another key ingredient in the operator design and implementation. 
A more detailed discussion of these issues can be found in Ref. [5]. 

A large effort was undertaken during the last two years to select optimal sets of baryon 
and meson operators in a large variety of isospin sectors and for zero momentum and 
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FIGURE 1. Hadron operator selection: low-statistics simulations have been performed to study the 
hundreds of single-hadron operators produced by our group-theoretical construction. A "pruning" pro- 
cedure was followed in each channel to select good sets of between six to a dozen operators. The plots 
above show the stationary-state energies extracted to date from correlation matrices of the finally selected 
single-hadron operators. Results were obtained using between 50 to 100 configurations on a 16^ X 128 
anisotropic lattice for Nf = 2 + 1 quark flavors with spacing ~ 0.12 fm, Us/at ~ 3.5, and quark masses 
such that niji ^ 380 MeV. Each box indicates the energy of one stationary state; the vertical height of each 
box indicates the statistical error 



non-zero on-axis, planar-diagonal, and cubic-diagonal momenta. Low-statistics Monte 
Carlo computations were done to accomplish these operator selections using between 
50 to 100 configurations on a 16^ x 128 anisotropic lattice for A^^ = 2 + 1 quark flavors 
with spacing as ~ 0. 12 fm, as /at ~ 3.5, and quark masses such that the pion has mass 
around 380 MeV. Stationary-state energies using the finally selected operator sets are 
shown in Fig. 1. The nucleon. A, S, E, and A baryons were studied, and light isovector 
and kaon mesons were investigated. Hundreds of operators were studied, and optimal 
sets containing eight or so operators in each symmetry channel were found. Future 
computations will focus solely on the operators in the optimal sets. 

A comprehensive picture of resonances requires that we go beyond a knowledge of 
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FIGURE 2. Diagrams of multi-hadron correlators that require having p noise sources on the later time 
t. Solution vectors are denoted by <p. (Left) A two-meson correlator. (Right) The correlator of a baryon- 
meson system. 



the ground state mass in each symmetry channel and obtain the masses of the low- 
est few states in each channel. This necessitates the use of matrices of correlation 
functions. Rather than evaluating a single correlator C{t), we determine a matrix of 
correlators Cij{t) = {Oi{t)0^j{tQ)) , where {OfJ = 1,...,A'^} are a basis of interpolat- 
ing operators with given quantum numbers. We then solve the generalized eigenvalue 
equation C{t)u = X{t,to)C{to)u to obtain a set of real (ordered) eigenvalues A„(f,fo), 
where Aq > Ai > ... > ^.n-i- At large Euclidean times, these eigenvalues then de- 
lineate between the different masses Xn{t,to) — > ^-^^"('-fo) ^ Q^^-AM„{t-to)-^^ where 
AM„ = min{| A/„ — Mi |: i ^ n}. The eigenvectors u are orthogonal with metric C(?o), 
and the eigenvectors yield information about the structure of the states. 

To study a particular eigenstate of interest with this method, all eigenstates lying 
below that state must first be extracted, and as the pion gets lighter in lattice QCD 
simulations, more and more multi-hadron states will lie below the excited resonances. 
A good baryon-meson operator of total zero momentum is typically a superposition of 
local interpolating fields at all sites on a time slice of the lattice. In the evaluation of the 
temporal correlations of such a multi-hadron operator, it is not possible to completely 
remove all summations over the spatial sites on the source time-slice using translation 
invariance. Hence, the need for estimates of the quark propagators from all spatial sites 
on a time slice to all spatial sites on another time slice cannot be sidestepped. Some 
correlators will involve diagrams with quark lines originating at the sink time t and 
terminating at the same sink time t (see Fig. 2), so quark propagators involving a large 
number of starting times t must also be handled. 

Finding better ways to stochastically estimate slice-to-slice quark propagators is cru- 
cial to the success of our excited-state hadron spectrum project at lighter pion masses. 
We have developed and tested a new scheme which combines a new way of smearing the 
quark field with a new way of introducing noise. The new quark-field smearing scheme, 
called Laplacian Heaviside (LapH), has been described in Ref. [4] and is defined by 

V/W = 0(ct2 + a)v/(x), (1) 

where A is the three-dimensional covariant Laplacian in terms of the stout-smeared 
gauge field and Os is the smearing parameter. The gauge-covariant Laplacian operator 
is ideal for smearing the quark field since it is one of the simplest operators that locally 



averages the field in such a way that all relevant symmetry transformation properties 
of the original field are preserved. Let V\ denote the unitary matrix whose columns 
are the eigenvectors of A, and let Aa denote a diagonal matrix whose elements are the 
eigenvalues of A such that A = VaA^V^. The LapH smearing matrix is then given by 

S = Va (of + Aa) V^. Let Vg denote the matrix whose columns are in one-to-one 
correspondence with the eigenvectors associated with the Nv lowest-lying eigenvalues 
of —A on each time slice. Then our LapH smearing matrix is well approximated by 
the Hermitian matrix S — Vg V/. Evaluating the temporal correlations of our hadron 
operators requires combining Dirac matrix elements associated with various quark lines 
Q. Since we construct our hadron operators out of covariantly-displaced, smeared quark 
fields, each and every quark line involves the following product of matrices: 

Q = D^j^SM-'^SD^''^\ (2) 

where D^'^ is a gauge-covariant displacement of type /. An exact treatment of such a 
quark line is very costly, so we resort to stochastic estimation. 

Random noise vectors rj whose expectations satisfy E{rii) = and ^(r/jr]*) = 5ij are 
useful for stochastically estimating the inverse of a large matrix M as follows[6]. Assume 
that for each of A';; noise vectors, we can solve the following linear system of equations: 
MZ('') = i]^'') for Then X^'^ = M~^ri^''\ and £:(X,T]T) = Mf.^ so that a Monte 

Carlo estimate of Mfj ^ is given by Mfj ^ fa lim^^^.c„ ^ £^ j X-^'^^ 7] j*^^ * . Unfortunately, this 
equation usually produces stochastic estimates with variances which are much too large 
to be useful. Variance reduction is done by diluting the noise vectors. A given dilution 
scheme can be viewed as the application of a complete set of projection operators P^'^^ . 

Define rj^^^ = ^k^'Oid: and further define zt"^! as the solution of Mi^xf^ = 77 H , then we 
have 

The use of Z4 noise ensures zero variance in the diagonal elements E{r]ir]*). 

The effectiveness of the variance reduction depends on the projectors chosen. With 
LapH smearing, noise vectors p can be introduced only in the LapH subspace. The noise 
vectors p now have spin, time, and Laplacian eigenmode number as their indices. Color 
and space indices get replaced by Laplacian eigenmode number. Again, each component 
of p is a random Z4 variable so that ^(p) = and ^(pp^) = /j. Dilution projectors P^*) 
are now matrices in the LapH subspace. In the stochastic LapH method, a quark line on 
a gauge configuration is estimated using 

where the subscripts m,v are compound indices combining space, time, color, and spin, 
and for a noise vector labelled by index r, displaced- smeared-diluted quark source and 
quark sink vectors can be defined by 



(5) 



(6) 



Our dilution projectors are products of time dilution, spin dilution, and Laph eigen- 
vector dilution projectors. For each type (time, spin, Laph eigenvector) of dilution, we 
studied four different dilution schemes. Let denote the dimension of the space of 
the dilution type of interest. For time dilution, N — Nt is the number of time slices on 
the lattice. For spin dilution, N = 4 is the number of Dirac spin components. For Laph 
eigenvector dilution, N = Ny is the number of eigenvectors retained. The four schemes 
we studied are defined below: 



where /, j = 0, . . . , — 1, and we assume N/K is an integer. We use a triplet (T, S, L) 
to specify a given dilution scheme, where "T" denote time, "S" denotes spin, and "L" 
denotes Laph eigenvector dilution. The schemes are denoted by 1 for no dilution, F for 
full dilution, and and for block-^ and interlace-^, respectively. For example, 
full time and spin dilution with interlace-8 Laph eigenvector dilution is denoted by 
(TF, SF, LIS). Introducing diluted noise in this way produces correlation functions with 
significantly reduced variances, yielding nearly an order of magnitude reduction in the 
statistical error over previous methods. The volume dependence of this new method was 
found to be very mild, allowing the method to be useful on large lattices. For all forward- 
time quark lines, we use dilution scheme (TF, SF, LIS), and for all same-sink-time quark 
lines, we use (Til 6, SF, LIS). 

Results for three isoscalar mesons are shown in Fig. 3. Such mesons are notoriously 
difficult to study in lattice QCD, but the new method appears to produce estimates 
of their temporal correlations with unprecedented accuracy. These plots suggest that 
evaluating correlation functions involving our multi-hadron operators will be feasible 
with the stochastic LapH method. 

We are currently carrying out these spectrum computations on 24^ x 12S and 32^ x 
256 anisotropic lattices with spatial spacing ~ 0.12 fm and aspect ratio as/ at ~ 3.5, 
where at is the temporal spacing, for pion masses mji ~ 380 MeV and m^r ~ 220 MeV. 
The calculations proceed in several steps: (a) generation of gauge-field configurations 
using the Monte Carlo method; (b) computation of quark sinks for various noises and 
dilution projectors using the configurations from the first step; (c) computation of the 
meson and baryon sources and sinks using the quark sinks from the second step; (d) 
evaluation of the correlators using the hadron sinks; (e) analysis of the correlators to 
extract the energies. Our results for the QCD stationary-state energies using, for the first 
time, both single-hadron and multi-hadron operators, should appear soon. 

This work was supported by the U.S. National Science Foundation under awards PHY- 
0510020, PHY-0653315, PHY-0704171, PHY-0969863, and PHY-0970137 and through 
TeraGrid resources provided by the Pittsburgh Supercomputer Center, the Texas Ad- 
vanced Computing Center, and the National Institute for Computational Sciences under 




5ij, a = 0, 

dij dau a = 0, 

^ij ^a, [Ki/N\ « = 0, 

^ij ^a,imoAK « = 0, 



, ^ — 1 , (interlace-^) 



,N-\ (full dilution) 
,K-\, (block-is:) 



(no dilution) 



80 



40 



-40 



1111111* 



iiii«« 



30 



20 



10 



CO 



60 
U40 
20 




10 20 



30 



10 20 
t 



30 



10 20 
t 



30 




5 10 15 20 5 10 15 20 

t t t 

FIGURE 3. Correlators C{t) and effective masses meff(f) against temporal separation t for single-site 
operators which produce the isoscalar pseudoscalar 77, vector o, and scalar a mesons. Results were 
obtained using 198 configurations with Nj — 2+\ flavors of quark loops on a 24^ x 128 anisotropic lattice 
with spacing ~ 0.12 fm and aspect ratio flj/a, ^ 3.5 for a pion mass iriji ~ 220 MeV. In the legends, 
"fwd" refers to contributions from the diagram containing only forward-time source-to-sink quark lines, 
"smt" refers to contributions from the diagram containing only quark lines that originate and terminate 
at the same time. For the <J channel, the "smt" contribution has a vacuum expectation value subtraction. 
Forward-time quark lines use dilution scheme (TF, SF, LI8) and same-time quark lines use (TI16, SF, 
LI8). 
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